{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import commpy as cp\n",
    "import scipy.signal as sig\n",
    "import scipy.linalg as la\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "from multiprocessing import Pool\n",
    "\n",
    "import traceback\n",
    "import warnings\n",
    "import sys\n",
    "\n",
    "%matplotlib inline\n",
    "DEFAULT_SEED = 100\n",
    "np.random.seed(DEFAULT_SEED)#set the random generator seed to default"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_topleitz(row_vector, channel_length):\n",
    "    first_row = row_vector[:channel_length]\n",
    "    first_col = row_vector\n",
    "    return la.toeplitz(first_col, first_row)\n",
    "\n",
    "def least_squares(A, b):\n",
    "    return np.linalg.lstsq(A, b)[0]\n",
    "\n",
    "def perform_least_squares(preamble, preamble_convolved, channel_length):\n",
    "    A, b = generate_topleitz(preamble, channel_length), preamble_convolved\n",
    "    prediction = least_squares(A, b)\n",
    "    return np.array(prediction)\n",
    "\n",
    "def add_awgn_noise(signal,SNR_dB):\n",
    "    \"\"\"  Adds AWGN noise vector to signal \n",
    "         to generate a resulting signal vector y of specified SNR in dB\n",
    "    \"\"\"\n",
    "    L=len(signal)\n",
    "    SNR = 10**(SNR_dB/10.0) #SNR to linear scale\n",
    "    Esym=np.sum(np.square(np.abs(signal)))/L #Calculate actual symbol energy\n",
    "    N0=Esym/SNR; #Find the noise spectral density\n",
    "    if(isinstance(signal[0], complex)):\n",
    "        noiseSigma=np.sqrt(N0/2.0)#Standard deviation for AWGN Noise when x is complex\n",
    "        n = noiseSigma*(np.random.randn(1,L)+1j*np.random.randn(1,L))#computed noise \n",
    "    else:\n",
    "        noiseSigma = np.sqrt(N0);#Standard deviation for AWGN Noise when x is real\n",
    "        n = noiseSigma*np.random.randn(1,L)#computed noise\n",
    "    y = signal + n #received signal\n",
    "    \n",
    "    return y.flatten()\n",
    "\n",
    "def modulate(data, mod_scheme='BPSK', demod=False):\n",
    "    \"\"\"  1. Modulates (or demodulates) data according to the modulation scheme \"\"\"\n",
    "    mod_schemes = ['BPSK', 'QPSK']\n",
    "    data = data.flatten()\n",
    "    if mod_scheme not in mod_schemes:\n",
    "        raise ValueError('Unknown modulation scheme, please choose from: '+ ' '.join(mod_schemes))\n",
    "    elif mod_scheme == 'QPSK':\n",
    "        modulator = cp.modulation.QAMModem(4)\n",
    "        if demod:\n",
    "            return modulator.demodulate(data, \"hard\")\n",
    "        return modulator.modulate(data)\n",
    "    elif mod_scheme == 'BPSK':\n",
    "        def bpsk_one(x):\n",
    "            if demod:\n",
    "                return 0 if x < 0 else 1\n",
    "            return -1 if x==0 else 1\n",
    "        bpsk = np.vectorize(bpsk_one)\n",
    "        return bpsk(data)\n",
    "    \n",
    "def apply_channel(signal, channel_function):\n",
    "    \"\"\"  2. Convolves signal with channel_function \"\"\"\n",
    "    channel_output = sig.convolve(signal, channel_function, mode='same') # convolve input complex data with the channel transfer function\n",
    "    return channel_output\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2000,) (2000,)\n",
      "[-1  1  1 ... -1  1  1]\n",
      "[-0.94827725  0.63083395  1.26572056 ... -0.63083395  0.63083395\n",
      "  1.26572056]\n",
      "[0.94827725 0.31744331]\n",
      "Predicted channel [0.94843661 0.31728395]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/ipykernel_launcher.py:7: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  import sys\n"
     ]
    }
   ],
   "source": [
    "channel_length = 2\n",
    "preamble_length = 2000\n",
    "data_length = 200\n",
    "\n",
    "modulation_scheme = 'BPSK'\n",
    "\n",
    "snr = 10\n",
    "\n",
    "channel_function = np.random.uniform(-1, 1, channel_length) \n",
    "channel_function = channel_function / np.linalg.norm(channel_function)\n",
    "\n",
    "# generate training data\n",
    "preamble_bits = np.random.randint(0,2, preamble_length) \n",
    "preamble_symbols = modulate(preamble_bits, modulation_scheme)\n",
    "preamble_convolved = apply_channel(preamble_symbols, channel_function)\n",
    "# generate testing data\n",
    "data_bits = np.random.randint(0,2, data_length)\n",
    "data_symbols = modulate(data_bits, modulation_scheme)\n",
    "data_convolved = add_awgn_noise(apply_channel(data_symbols, channel_function), snr)\n",
    "print(preamble_symbols.shape, preamble_convolved.shape)\n",
    "print(preamble_symbols)\n",
    "print(preamble_convolved)\n",
    "print(channel_function)\n",
    "channel_estimate = perform_least_squares(preamble_symbols, preamble_convolved, channel_length)\n",
    "print('Predicted channel', channel_estimate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_A_hat(channel_estimate, data_length):\n",
    "    matrix = np.zeros((data_length, data_length))\n",
    "    for i in range(data_length - 1): \n",
    "        matrix[i][i] = channel_estimate[0]\n",
    "        matrix[i+1][i] = channel_estimate[1]\n",
    "    matrix[data_length - 1][data_length - 1] = channel_estimate[0]\n",
    "    return matrix\n",
    "\n",
    "def calc_accuracy(original, predictions):\n",
    "    return np.sum(original == predictions) / (original.shape[0])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.94843661 0.         0.         ... 0.         0.         0.        ]\n",
      " [0.31728395 0.94843661 0.         ... 0.         0.         0.        ]\n",
      " [0.         0.31728395 0.94843661 ... 0.         0.         0.        ]\n",
      " ...\n",
      " [0.         0.         0.         ... 0.94843661 0.         0.        ]\n",
      " [0.         0.         0.         ... 0.31728395 0.94843661 0.        ]\n",
      " [0.         0.         0.         ... 0.         0.31728395 0.94843661]]\n",
      "RHS\n",
      "[[ 0.94843661 -0.63115266  0.63115266 ...  1.26572056  1.26572056\n",
      "   1.26572056]\n",
      " [-0.94843661  0.63115266 -0.63115266 ... -1.26572056 -1.26572056\n",
      "  -1.26572056]\n",
      " [ 0.94843661 -0.63115266  0.63115266 ...  1.26572056  1.26572056\n",
      "   1.26572056]\n",
      " ...\n",
      " [ 0.94843661 -0.63115266  0.63115266 ...  1.26572056  1.26572056\n",
      "   1.26572056]\n",
      " [ 0.94843661 -0.63115266  0.63115266 ...  1.26572056  1.26572056\n",
      "   1.26572056]\n",
      " [ 0.94843661 -0.63115266  0.63115266 ...  1.26572056  1.26572056\n",
      "   1.26572056]]\n",
      "Inverse\n",
      "[[ 9.9543455   0.03038153 -0.03038153 ... -0.06092746 -0.06092746\n",
      "  -0.06092746]\n",
      " [ 0.03038153  9.97978211  0.02021789 ...  0.04054518  0.04054518\n",
      "   0.04054518]\n",
      " [-0.03038153  0.02021789  9.97978211 ... -0.04054518 -0.04054518\n",
      "  -0.04054518]\n",
      " ...\n",
      " [-0.06092746  0.04054518 -0.04054518 ...  9.91869025 -0.08130975\n",
      "  -0.08130975]\n",
      " [-0.06092746  0.04054518 -0.04054518 ... -0.08130975  9.91869025\n",
      "  -0.08130975]\n",
      " [-0.06092746  0.04054518 -0.04054518 ... -0.08130975 -0.08130975\n",
      "   9.91869025]]\n",
      "Right inverse\n",
      "[[ 9.9543455   0.03038153 -0.03038153 ... -0.06092746 -0.06092746\n",
      "  -0.06092746]\n",
      " [ 0.03038153  9.97978211  0.02021789 ...  0.04054518  0.04054518\n",
      "   0.04054518]\n",
      " [-0.03038153  0.02021789  9.97978211 ... -0.04054518 -0.04054518\n",
      "  -0.04054518]\n",
      " ...\n",
      " [-0.06092746  0.04054518 -0.04054518 ...  9.91869025 -0.08130975\n",
      "  -0.08130975]\n",
      " [-0.06092746  0.04054518 -0.04054518 ... -0.08130975  9.91869025\n",
      "  -0.08130975]\n",
      " [-0.06092746  0.04054518 -0.04054518 ... -0.08130975 -0.08130975\n",
      "   9.91869025]]\n",
      "Data estimate\n",
      "[ 1.00756829 -1.00756829  1.00756829 -1.00756829  1.00756829  1.00756829\n",
      " -1.00756829 -1.00756829 -1.00756829 -1.00756829  1.00756829  1.00756829\n",
      " -1.00756829  1.00756829 -1.00756829 -1.00756829 -1.00756829  1.00756829\n",
      " -1.00756829  1.00756829  1.00756829 -1.00756829 -1.00756829 -1.00756829\n",
      "  1.00756829 -1.00756829  1.00756829 -1.00756829  1.00756829 -1.00756829\n",
      "  1.00756829 -1.00756829 -1.00756829 -1.00756829  1.00756829 -1.00756829\n",
      "  1.00756829  1.00756829 -1.00756829 -1.00756829  1.00756829  1.00756829\n",
      " -1.00756829 -1.00756829  1.00756829  1.00756829 -1.00756829 -1.00756829\n",
      " -1.00756829 -1.00756829 -1.00756829 -1.00756829  1.00756829 -1.00756829\n",
      " -1.00756829  1.00756829  1.00756829  1.00756829  1.00756829  1.00756829\n",
      "  1.00756829  1.00756829 -1.00756829 -1.00756829  1.00756829 -1.00756829\n",
      "  1.00756829 -1.00756829 -1.00756829 -1.00756829 -1.00756829 -1.00756829\n",
      "  1.00756829 -1.00756829  1.00756829  1.00756829  1.00756829 -1.00756829\n",
      "  1.00756829 -1.00756829  1.00756829 -1.00756829 -1.00756829  1.00756829\n",
      " -1.00756829  1.00756829 -1.00756829  1.00756829 -1.00756829  1.00756829\n",
      "  1.00756829 -1.00756829  1.00756829  1.00756829 -1.00756829  1.00756829\n",
      "  1.00756829 -1.00756829  1.00756829 -1.00756829 -1.00756829 -1.00756829\n",
      " -1.00756829 -1.00756829  1.00756829  1.00756829 -1.00756829 -1.00756829\n",
      " -1.00756829  1.00756829 -1.00756829 -1.00756829 -1.00756829  1.00756829\n",
      "  1.00756829 -1.00756829 -1.00756829 -1.00756829  1.00756829  1.00756829\n",
      "  1.00756829 -1.00756829  1.00756829  1.00756829  1.00756829 -1.00756829\n",
      "  1.00756829  1.00756829 -1.00756829 -1.00756829 -1.00756829  1.00756829\n",
      "  1.00756829 -1.00756829 -1.00756829 -1.00756829 -1.00756829 -1.00756829\n",
      " -1.00756829  1.00756829 -1.00756829  1.00756829 -1.00756829  1.00756829\n",
      " -1.00756829 -1.00756829  1.00756829 -1.00756829 -1.00756829 -1.00756829\n",
      "  1.00756829 -1.00756829  1.00756829 -1.00756829 -1.00756829  1.00756829\n",
      "  1.00756829  1.00756829 -1.00756829 -1.00756829 -1.00756829  1.00756829\n",
      " -1.00756829  1.00756829  1.00756829  1.00756829  1.00756829  1.00756829\n",
      "  1.00756829  1.00756829  1.00756829  1.00756829  1.00756829 -1.00756829\n",
      " -1.00756829 -1.00756829 -1.00756829 -1.00756829 -1.00756829 -1.00756829\n",
      " -1.00756829 -1.00756829  1.00756829 -1.00756829  1.00756829 -1.00756829\n",
      " -1.00756829  1.00756829  1.00756829  1.00756829 -1.00756829 -1.00756829\n",
      "  1.00756829 -1.00756829  1.00756829 -1.00756829  1.00756829  1.00756829\n",
      "  1.00756829  1.00756829]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 128,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_hat = generate_A_hat(channel_estimate, data_length)\n",
    "print(A_hat)\n",
    "# print('Ahat')\n",
    "# print(A_hat)\n",
    "# print('data convolved')\n",
    "# print(data_convolved)\n",
    "# print('data convolved 2')\n",
    "# data_convolved_estimate = np.dot(A_hat, data_symbols)\n",
    "# print(data_convolved_estimate)\n",
    "# print('Norm', np.linalg.norm(data_convolved_estimate - data_convolved))\n",
    "std = 10 ** (-snr / 20)\n",
    "outer_prod_data = np.outer(data_symbols, data_symbols) #xxT\n",
    "RHS = np.dot(outer_prod_data, A_hat.T) # xx_TA_T\n",
    "to_invert = np.dot(np.dot(A_hat, outer_prod_data), A_hat.T) + (std ** 2) * np.identity(A_hat.shape[0]) #A_hat x x.T A_hatT\n",
    "B = np.dot(RHS, np.linalg.inv(to_invert))\n",
    "print('RHS')\n",
    "print(RHS)\n",
    "print('Inverse')\n",
    "print(np.linalg.inv(to_invert))\n",
    "print('Right inverse')\n",
    "print(right_inv)\n",
    "# B = np.dot(RHS, np.identity(data_length))\n",
    "data_est = np.dot(B, data_convolved)\n",
    "data_symbols_est = np.where(data_est > 0, 1, -1)\n",
    "print('Data estimate')\n",
    "print(data_est)\n",
    "calc_accuracy(data_symbols_est, data_symbols)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(200, 200)\n",
      "(200,)\n"
     ]
    }
   ],
   "source": [
    "print(A_hat.shape)\n",
    "print(data_symbols.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
